On the lack of universality in decaying magnetohydrodynamic turbulence 
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Using computations of three-dimensional magnetohydrodynamic (MHD) turbulence with a 
Taylor-Green flow, whose inherent time-independent symmetries are implemented numerically, and 
in the absence of either a forcing function or an imposed uniform magnetic field, we show that three 
different inertial ranges for the energy spectrum may emerge for three different initial magnetic 
fields, the selecting parameter being the ratio of the Alfven to the eddy turnover time. Equivalent 
computational grids range from 128 3 to 2048 3 points with a unit magnetic Prandtl number and a 
Taylor Reynolds number of up to 1500 at the peak of dissipation. We also show convergence of our 
results with Reynolds number. Our study is consistent with previous findings of a variety of energy 
spectra in MHD turbulence by studies performed in the presence of both a forcing term with a given 
correlation time and a strong, uniform magnetic field. In contrast to the previous studies, however, 
the ratio of characteristic time scales here can only be ascribed to the intrinsic nonlinear dynamics 
of the flows under study. 

PACS numbers: 47.65.-d; 47.27.Jv; 94.05.Lk; 95.30.Qd 



I. INTRODUCTION 

Turbulence forms the backbone of many natural phe- 
nomena in the atmosphere and ocean, as well as in as- 
trophysical flows. In the latter, it is often accompanied 
by the coupling of vortices and current structures. For 
incompressible neutral fluids, under the assumption of 
a high Reynolds number (and therefore a long dissipa- 
tion time), as is the case for many geophysical flows, 
the relevant time scale for the problem is the nonlinear 
eddy turnover time. For such flows, the phenomenol- 
ogy developed by Kolmogorov in 1941 (hereafter referred 
to as "K41") [1(, predicting a kinetic energy spectrum 
E^ A1 {k) ~ fc -5 / 3 , represents a good first approach even if 
corrections to this phenomenology for higher-order statis- 
tics are known to exist, due to the breakdown of the 
self-similarity represented by a simple power-law energy 
spectrum. When electromagnetic forces are introduced, 
other time scales can arise, such as the Alfven time, as- 
sociated with the propagation of transverse waves along 
magnetic field lines. K41 phenomenology may still ap- 
ply, but one must also consider the role of Alfven waves 
in producing a different power law for the total energy 
spectrum, as illustrated independently by Iroshnikov and 
Kraichnan (hereafter, "IK") 0: Ef K (k) ~ fc~ 3 / 2 . 

Isotropy is assumed by both K41 and IK, but it is not 
necessarily achieved. In neutral flows, if the anisotropy 
of the small scales, in the form of elongated vortex fila- 
ments, occurs locally in space, isotropy may be recovered 
overall because the filaments are randomly oriented, and 
the vorticity spectrum k 2 E v (fc), which peaks in the small 
scales, contributes little to the large-scale spectrum. In 
contrast, the anisotropy of MHD originates from a large- 



scale magnetic field, which can be dominant energetically 
and relevant at all scales. Studies of anisotropic MHD 
date back to the mid-fifties for liquid metals at low mag- 
netic Reynolds number Q (see alsofi |) , and a bit later for 
fully developed MHD turbulence [a, la] • More recently, a 
wealth of new studies on MHD turbulence has been made 
possible @, HJjl in part by the revival of weak turbulence 
theory (e.g. [H| for MHD), the availability of more de- 



numerical simulations [1 



tailed observations [111, 12|, and improved resolution in 
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From the theoretical point of view, the presence of 
a strong background magnetic field Bo allows for the 
existence of a small parameter characterizing the ratio 
of velocity and magnetic field fluctuations to |Bo|, en- 
abling an analytical solution via the weak turbulence 
(hereafter, "WT") approach. In contrast to the K41 
and IK spectra, the WT energy spectrum is anisotropic: 
EwT(kj_, fc|| ) ~ kj 2 , where perpendicular and parallel 
are relative to the direction of Bo; there is, in fact, no 
prescribed transfer in fen at lowest order. It is interesting 
to note that the IK phenomenology is compatible with 
weak turbulence in the isotropic limit k± re fcy, giving it 
a stronger theoretical footing. Other phenomenological 
approaches hypothesize that even with a strong back- 
ground field Bo making the flow highly anisotropic, an 
"anisotropic Kolmogorov" scaling of the energy spectrum 
is possible by way of a dynamical effect that makes the 
two characteristic times of the problem (the Alfven time 
and the eddy turnover time) equal at all scales [l8j|. In 
fact, when this hypothesis is relaxed to a constant ratio 
across the inertial range (not necessarily equal to unity) , 
the dynamics can then be shown to be compatible with a 
variety of inertial range scalings, including K41, IK and 
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WT [8| . Weak MHD turbulence has been observed in the 
magnetosphere of Jupiter [Til ] . where the strong Jovian 
field creates a favorable environment for wave interac- 
tions to dominate the dynamics. In the solar wind, data 
for a long time have indicated that the spectrum appears 
Kolmogorovian [l9j |. although recent observations indi- 
cate a more complex dynamics (discussed further in the 
conclusion, fjV|. 

Numerical simulations to date are unable to give a 
definitive answer to the question of spectral index in 
plasma turbulence, at least in three dimensions. The dif- 
ference between the K41 and IK scalings is subtle enough 
that any type of contamination, in particular by intermit- 
tency as well as dissipative small-scale effects, will blur 
the results. Intcrmittency, the sporadic occurrence of in- 
tense small-scale structures, tends to steepen the direct 
cascade energy spectrum. In fact, it has been shown 
both in two and three spatial dimensions that intermit- 
tency in MHD generally leads to stronger corrections 
to high-order structure functions than in neutral fluids 
[3 HOj HH , and that the magnetic field is more intermit- 
tent than the velocity [2(| H2] • 

Recent studies indicate that, in the presence of an ex- 
ternal force with a given autocorrelation time and/or a 
strong, uniform magnetic field, the ener gy s pectrum can 
exhibit different power laws (e.g., [TEl. l23l. [24[). Such 
varied spectral indices can be ascribed to the variation 
of time scales [l3j] or to the presence of complex struc- 
tures, such as ribbons (see [151 ] and references therein, 
[25l]). However, the possibly simpler problem of incom- 
pressible MHD decay in three dimensions with Bo = 
has not been examined in this light [2(|. Therefore, it 
is the purpose of this paper to do so by way of high- 
resolution numerical simulation and to show that indeed 
several classes of dynamics are possible in decaying MHD 
turbulence. In the next section, we give equations and 
initial conditions; EIIIII is dedicated to the temporal be- 
havior of the flows, mV\ to the spectra observed in this 
paper and convergence of the results with Reynolds num- 
ber, and |jV]to a discussion and brief concluding remarks. 



II. THREE CLASSES OF MAGNETIC 
TAYLOR-GREEN FLOW 

The MHD equations for an incompressible flow with a 
velocity v and magnetic induction b (in units of Alfvcn 
velocity) read: 



total energy E T = (v 2 +b 2 )/2 = E V +E M , cross helicity 



dv 



V'Vv = -Vp + jxb + i/Av, (1) 



db 

~dt 



V x (v x b) + ^Ab , 



(2) 



V- v = = V-b , 



with p the fluid pressure and j = V x b the current 
density. In the absence of viscosity v and resistivity r], the 



(v • b}/2, and magnetic helicity H 



M 



b)/2 



(where b = V x a defines the magnetic potential, a) 
are conserved. The Reynolds number here is defined as 
Re = v rms L T /is, and the magnetic Reynolds number as 
Rm = v rms L M /r], with the integral scale and kinetic and 
magnetic integral scales defined, respectively, as: 



L V,M,T = ^ 



J E V ' M > T (k)k~ 1 dk 
jE V ' M < T (k)dk 



Similarly, one can define the Taylor Reynolds number 
R\ = v rms \ T /i/, where the Taylor scale A T is defined as 



X 1 =2tt 



fE T {k)dk 



1/2 



fE T (k)k 2 dk, 

Kinetic and magnetic Taylor scales can also be defined: 



\ V < M = 27T 



fE v > M (k)dk 
fE v > M {k)k 2 dk 



1/2 



Ideal MHD {u = 0, T) = 0) has been studied numeri- 
cally both in two dimensions [28| and in three [29L [30I ] , 
including with adaptive mesh refinement [3l| . Such sim- 
ulations are important for understanding the initial non- 
linear exchanges among modes until the smallest resolved 
scales are reached, at which point dissipation must be 
introduced to continue the computation and reach a 
fully developed turbulent flow with current and vortic- 
ity sheets. 

The velocity field we choose for our initial conditions 
is the Taylor-Green (hereafter, "TG") vortex (3^1 corre- 
sponding to a von Karman flow between two counter- 
rotating disks. The simplest TG velocity field can be 
written as (33[ (see also [HI): 

, WTG(x 1 y, z) = vq [sin x cos y cos z e x — cos x sin y cos z e y ] 

The velocity component in the third direction is zero 
initially but grows with time. We also define ujtg = 
V x vtg (and u> = V x v, the vorticity as usual). This 
TG vortex has been used not only in numerical studies 
[35| . but also extensively in laboratory studies of fluid 
turbulence and as a driver for the generation of magnetic 
fields within liquid metals flows [361 ]. 

A generalization of the TG vortex symmetries to MHD 
was presented in [3fj| . where the ideal case was studied 
with the following initial magnetic field configuration: 



bl = 

bl = 



60 cos(x) sin(y) sin(z) 
6q sin(x) cos(y) sin(z) 
— 2&Q sin(x) sin(y) cos(z). 



(3) 
(4) 

(5) 



It can be shown that the magnetic field h 1 is every- 
where perpendicular to the faces, or "walls," of a box de- 
fined as [0, 7r] 3 . The current j 7 = V x h 1 is then found to 
be parallel to the walls, which thus can be considered as 
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electrical insulators; for this reason, we refer to this type 
of TG flow as "insulating." Also of interest is that the 
cross-correlation He between the velocity and magnetic 
fields is identically zero globally and h 1 = — \bo/vo\u>TG- 
Another magnetic field was proposed in [301 ] . namely: 

b° = = b% sm(2x) cos(2y) cos(2z) (6) 
b c y = = b% cos(2.t) sin(2y) cos(2z) (7) 
b C z = = -2&£cos(2x)cos(2y)sin(2z) , (8) 

where the current in the [0, 7r] 3 box here is perpendicular 
to the walls, which are therefore "conducting." In this 
configuration, H c is non-zero but weak (less than 4% 
at its maximum over time, in a dimcnsionlcss measure 
relative to the total energy). 

Finally, we also introduce an alternative to the insu- 
lating magnetic induction above, which we call h A (for 
"alternative" insulating flow), defined as: 

b£ = = &o cos(2z) sin(2y) sin(2^) (9) 
by = = -&o sin(2a;) cos(2y) sin(2z) (10) 
bf = =0, (11) 

configuration for which again He = 0. Note that the 
magnetic helicity is zero in all three configurations. 

All flows are initialized at the largest scales in order 
to obtain the highest possible Reynolds number, and all 
have unit magnetic Prandtl number (i.e., v = rj). The 
values of the parameters for all runs described in this 
paper are given in Table fl] For the three classes of flow 
proposed here (named hereafter "I," "A," and "C," re- 
ferring to the three induction configurations), the same 
TG velocity is applied at t = 0, and the fields are nor- 
malized such that E v (t = 0) = E M (t = 0) = 0.125: at 
the start of each simulation, the kinetic and magnetic 
energies are in equipartition. The resolutions of the runs 
range from 64 3 to 2048 3 points (in terms of equivalent 
grids for computations that would not exploit the sym- 
metries; see [30(1), allowing the range of Reynolds number 
to vary over a factor close to 40. 

Taylor-Green configurations possess several inherent 
symmetries within a cube of length 27r (periodicity). Mir- 
ror symmetries about the planes x = 0, x = tt, y = 0, 
y = 7r, z = 0, and z = it; and rotational symmetries of an- 
gle Ntt about the axes (x,y, z) = (§,y, f) and (x, §,§) 
and of angle Ntt/2 about the axis (5, ^, z), for N G Z. 
All flows are defined in the [0, 27r] 3 box and satisfy pe- 
riodic boundary conditions of the domain. Within the 
domain, the planes of mirror symmetry mentioned above 
form the insulating and conducting walls of the smaller 
[0,7r] 3 boxes. 

As the symmetries of the TG flows are also symmetries 
of the MHD equations, they are preserved by time evo- 
lution of the solutions. Numerical implementation of the 
symmetries allows for substantial savings in both com- 
puting time and memory usage at a given Reynolds num- 
ber, with no approximation or closure scheme needed. 



The numerical method is pseudo-spectral, with mini- 
mum wavenumber k m i n = 1 and maximum wavenumbcr 
kmax = N/3, where N the number of grid points in each 
direction, using a 2/3-deliasing rule; the temporal inte- 
gration is performed using a second-order Runge-Kutta 
scheme. The code follows the parallelization, using MPI, 
as described in [33] (see [33| for a detailed implementa- 
tion for Euler and Navier-Stokcs and [3(| for MHD). It 
was checked for runs at a resolution of 512 3 grid points 
that the differences between results obtained with the 
code implementing all the symmetries and those with a 
general code (in which the symmetries are not imposed 
explicitly, but the initial conditions are the same) were 
small and did not grow in time: in Fig. [T] is shown the 
relative difference for the domain-integrated total enstro- 
phy, defined as (uj 2 + j 2 ) /2, which remains of the order 
of 10~ 5 throughout the computation (see also [3fJ). Fur- 
ther results concerning such a comparison, as well as an 
analysis of the overall dynamics of the 16 flow are given 




2 4 6 8 

t 



FIG. 1: Relative error in the total enstrophy as a function 
of time for the 13 flow (Re = 940) between a full direct nu- 
merical simulation on a grid of 512 3 points with no imposed 
symmetries and one making use of code-implemented symme- 
tries using the same grid. 



III. THE EMERGENCE OF DIFFERENT 
REGIMES 

In a three-dimensional non-rotating turbulent flow, 
the nonlinear terms — here, the momentum advection, 
Lorentz force, and induction terms — produce coupling 
between modes, and both the kinetic and magnetic en- 
ergy are transferred to small scales as a result. As pro- 
gressively smaller scales become excited, the volume- 
integrated vorticity and current density grow until dis- 
sipation sets in, as can be seen in Fig. [5] (top) for the 
highest-resolution runs 16 (solid line), A6 (dashed line) 
and C6 (dotted line), following the nomenclature of Ta- 
ble fl] The total level of dissipation is comparable for 
runs A6 and C6 but substantially lower for run 16. The 
ratios of magnetic to kinetic dissipation are smaller than 
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FIG. 2: ( Top) Temporal evolution of the total dissipation 
for the highest Reynolds numbers for each type of flow, rep- 
resented by runs 16 (solid), A6 (dashed), and C6 (dotted), 
as described in Table HI (Bottom) Ratio of total magnetic to 
kinetic energy E /E for these same runs. Note that 16 has 
noticeably less dissipation and more magnetic energy. 



their equivalent energy counterparts, and they reach their 
maxima earlier (not shown) ; this indicates a stronger ten- 
dency towards cquipartition in the small scales, as can be 
expected since the Alfven time varies as 1/k (the data for 
run 16 is discussed in [l7j|). 

The energy exchanges between the velocity and mag- 
netic fields are complex, yielding different ratios of mag- 
netic to kinetic energy that change both in time and from 
flow to flow, as we see in Fig. [2] (bottom). These differ- 
ences can also be understood in terms of the diversity 
of nonlinear terms in the MHD equations, their relative 
importance depending on the initial conditions. For ex- 
ample, runs I arc dominated by magnetic energy after 
a short period of time, although the initial fields are in 
equipartition. For this flow, the magnetic field and the 
vorticity are initially parallel at every point in space. As 
a result, the nonlinear terms in the MHD equations ini- 
tially lead to a more rapid production of magnetic energy 
than in the other flows. At later times, the action of the 
Lorentz force differentiates the evolution of the magnetic 
field from that of the vorticity. 

We find thus that the intrinsic nonlinear dynamics for 
the three sets of initial conditions lead to magnetically 
dominated (I) flows, quasi-cquipartitioned (A) flows, or 
flows with sub-dominant magnetic energy (C), at least 



TABLE I: Parameters used in the simulations. N is the equiv- 
alent linear grid resolution; v rms and b rma are the rms veloc- 
ity and magnetic field at peak of dissipation; v = r\ is the 
kinematic viscosity; and Re is the Reynolds number based on 
the integral scale of the flow. Note the growth of brms with 
Re in all but one example, as well as the diminishing ratio 
brms/ Vrms going from I to A to C flows. * 13 has also been 
run on a full grid, for comparison purposes (see Fig. [TJ. 
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FIG. 3: Ratio of magnetic to kinetic energy spectra averaged 
over At = 0.5 (1.5 to 2 turnover times) about the maximum of 
dissipation for run 16 (solid), A6 (dashed), and C6 (dotted). 
(Labels are the same as in Fig. [2] see Table U for run details.) 
Note the dominance of the magnetic energy at large scale for 
run 16 (solid) , and the tendency towards equipartition at small 
scales for all runs, with a slight excess of magnetic energy. 

in the large scales. Fig. [3] displays the variation with 
wavenumber of the ratio of the magnetic to kinetic energy 
spectra for the three high-resolution runs, one for each 
class. A surplus of magnetic energy at large scale for the 
1 flow is evident, as is a tendency in all flows towards 
quasi-cquipartition at small scales. A slight excess of 
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FIG. 4: Total energy spectra (top) compensated by fc 5//3 and 
averaged over At = 0.5 (1.5 to 2 turnover times) about the 
maximum of dissipation, and ratio of nonlinear to Alfven time 
scales as a function of wavenumber (bottom) for the same runs 
(labels as in Fig. [2] solid line for I, dash for A and dots for C). 
Slopes are given only as a reference. The three arrows indicate 
the magnetic Taylor scale. Note that the three spectra follow 
noticeably different spectral laws and possibly different scale- 
dependence for their time scales as well (see text). 



magnetic excitation at small scales can be observed, a 
feature perhaps linked to the absence of a magnetically 
induced "eddy viscosity" for the magnetic energy akin to 
one for the kinetic energy, as shown in studies of MHD 
turbulence using transport coefficients derived from the 
full intcgro-differential equations arising from a two-point 
closure [381 ]. 

When examining now the resulting energy spectra, av- 
eraged over an interval of time At = 0.5 about the 
maximum of dissipation of each flow, one can easily 
distinguish the three flows, with measurably different 
power laws. Fig. [J] (top) gives the total energy spec- 
tra for the same three runs, compensated by fc 5 / 3 , cal- 
culated from data averaged over eleven temporal out- 
puts at t G [3.75, 4.25] for 16, t £ [4.5, 5.0] for A6 and 
t E [4.75, 5.25] for C6. The A6 flow is near the K41 
scaling; the C6 flow has a shallower spectrum close to 
the IK dynamics, with a fc~ 3 / 2 index; and the 16 flow has 
a steeper spectrum close to the WT expectation, with a 
k~ 2 power law (see also Section HVl ). The denotations of 
the spectra as K41, IK, and WT are used here for sim- 
plicity, and as will be discussed later, more simulations 



arc needed to decide whether these arc the dynamical at- 
tractors of the equations or if other solutions exist. Re- 
gardless, the three sets of initial conditions clearly lead 
to different spectral behavior, which can be linked to the 
several times scales involved in the system as shown next. 
We recall that in [l4[ , the IK spectrum was followed for 
larger wavenumbcrs by a steeper spectrum, kj 2 corre- 
sponding to WT. The data in the present study are not 
quite sufficient to confirm this finding, but it is still highly 
possible that the IK spectrum in C6 (dotted line in Fig. 2]) 
is followed by a steeper behavior, with a transition occur- 
ring, as in |14| . at the magnetic Taylor scale, indicated 
by an arrow, whereas no such transition is visible for the 
other two classes, A and I. 

In Fig. [5] (bottom) is shown the ratio of the nonlinear 
time to the Alfven time, R T (k) = r e (fc)/r J 4(fc), where 
these times are defined as: 

re(fc) = -^wwm ' TA[k) = ik ' 

with Bo the field in the largest-scale mode; hence: 

Rr{k) = B » . 

v / 2kE v (k) 

Contrasting the plots in Fig. [4l we may conjecture 
that it is the competition between the two characteristic 
phenomena in MHD turbulence (nonlinear steepening as 
measured by r e , and wave interactions as measured by 
ta) that produces different equilibria among scales and 
therefore different energy spectra. The fact that these ra- 
tios here are always greater than unity is not significant 
in itself, as phenomenological determination of charac- 
teristic times leaves them within some constant factor of 
order unity, but what may be significant is their variation 
from flow to flow, as well as their variation with scale. 

As discussed earlier, it has been hypothesized that 
MHD turbulence dynamics may be understood in the 
context of an equilibrium between turbulent eddies and 
Alfven wave interactions @, H, H, [l|[ . Indeed, the nonlin- 
ear MHD equations accept the solutions v = ±b. MHD 
turbulence can also be viewed as the competition between 
nonlinear steepening and (semi)-dispersive effects, some- 
what akin to soliton dynamics [23|. The nonlinear com- 
petition between eddies and waves, in this light, perhaps 
could be measured by R T . By setting R T = 1 at all scales, 
a K41-type spectrum could occur; relaxing the condition 
and leaving R T equal to a constant, independent of scale, 
a model can be devised [l8j whereby such a condition can 
be made compatible with the IK phenomenology and WT 
theory as well (Hi- 
lt could be argued that the hypothesis of constancy of 
R T with scale is verified by the data displayed in Fig. [4] 
However, if we estimate as usual the nonlinear time as 
l/- s /2k 3 E v (k), then R T must vary as fc 1 / 3 , fc 1 / 4 , and 
fc 1 / 2 , respectively, for K41, IK, and WT dynamics. The 
results shown in Fig. [5] are also compatible with this in- 
terpretation, favoring a slight steepening of R T (k) as we 
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go from I to A to C initial conditions. This observa- 
tion may be affected, though, by a contamination of the 
inertial range by dissipation associated with still some- 
what moderate Reynolds numbers. Obviously, higher 
Reynolds number computations will help to clarify this 
point. 
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FIG. 5: Variation with Reynolds number Re of the ratio of 
eddy turn-over to Alfven time scales computed at the Taylor 
scale A A/ for each flow. 11-16 (+), A1-A6 (o), and C1-C6 (*) 
are plotted in order of resolution and Reynolds number, as 
listed in Table [T] All quantities are computed in an interval of 
At = 0.5 about the peak of dissipation for each run. Note the 
rather different values of these ratios and onset of saturation 
for the highest Re indicative of the beginning of a convergence 
to a high-f?e regime. 



TABLE II: Integral scales L T ' V ' M , Taylor scales \ V ' M , and 
Taylor Reynolds number R^ based on A T ; all values are taken 
at the instantaneous peak of dissipation. Note the substan- 
tially larger integral and Taylor scales in run 16, as well as the 
kinetic integral scale L v in C6. These scales directly affect 
the Reynolds numbers of the flows and the relative time scale 
responsible for the dynamics. 
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IV. SCALING WITH REYNOLDS NUMBER 

It is notoriously difficult to measure spectral indices 
of power laws found in numerical simulations, in partic- 
ular because of the small extent of the inertial range, 
sandwiched in wavenumber space between the energy- 
containing and dissipative scales. A further difficulty 
arises when more than one inertial regime exists, as 
found, for example, in Hall MHD [3i| and in weak turbu- 
lence [14| . The question then arises as to whether or not 
the spectra plotted in Fig. 0] and other data presented 
here would be consistent with the dynamics of equiva- 
lent flows at substantially higher Reynolds numbers. To 



address such issues we now turn to a convergence study 
of the data. 

In Fig.[U R T is plotted as a function of scale, but to see 
its dependence on Reynolds number, we can compute it 
for a specific scale in the inertial range, which we choose 
to be the Taylor scale A M . Fig. shows how R T {\ M ) 
changes with Reynolds number for each class of initial 
conditions. We see that this ratio tends towards a con- 
stant whose value depends on the type of flow, although 
higher Reynolds numbers should be investigated to con- 
firm this tendency. The transition to a converged value 
of R T appears to take place at Re « 3 x 10 3 , after a 
long evolution over a sequence of smaller Reynolds num- 
bers. For this reason, simulations at moderate Reynolds 
numbers can only tentatively show the asymptotic results 
needed to understand fully developed turbulence. 

Other indicators can be examined to further suggest 
convergence, at least for the highest two resolutions (runs 
5 and 6) for each class of flow. In Fig. [()] we show — again 
as a function of Reynolds number — the variation of the 
maximum of dissipation (top), the time at which this 
maximum is reached (middle), and the Taylor Reynolds 
number (bottom). We observe that the maximum of 
dissipation for the A and C initial conditions tends to 
level off towards a constant value as the Reynolds num- 
ber is increased, as seen before in two-dimensional MHD 
pfl . l4l| . three-dimensional Navier-Stokes [12], and an 
earlier three-dimensional MHD study (43|. The I flows 
seem to show a different trend, so higher resolution runs 
of this class (at higher Reynolds numbers) are necessary 
to determine if the maximum of dissipation docs indeed 
reach a constant; here we can only observe a slowing of its 
decrease, and in fact the trend follows a power-law when 
examining the first maximum (as opposed to the absolute 
maximum, at f ~ 3.25 and « 4.75 respectively, see Fig. 
[2 top and [13). It should nevertheless be noted that, 
for each class, the maximum of dissipation occurs at a 
time that depends significantly on the Reynolds number 
(Fig. [TjJ middle). Finally, we observe a clear scaling of 
the Taylor Reynolds number R\ [44[ with the large-scale 
Reynolds number Re (bottom): the two flows at the high- 
est Reynolds numbers for each class are consistent with 
Rx ~ Re 1 / 2 , indicating again that we have reached an 
asymptotic state. 

With reasonable evidence of convergence of flows be- 
yond a Reynolds number threshold, we can now turn to 
the energy spectra we observe for each class of flow. Fig. [7] 
displays, for each class of initial conditions, the total en- 
ergy spectra averaged again over an interval At = 0.5 
about the maximum of dissipation of each run (see figure 
caption and Table UJ. The spectra arc compensated by 
k 2 for the I flows, by k 5 / 3 for the A flows, and by fc 3 / 2 for 
the C flows. These plots strongly suggest that the scaling 
predictions of WT, K41, and IK give credible descrip- 
tions of the I (top), A (middle), and C (bottom) flows, 
respectively, recognizing that intermittency can steepen 
the spectrum of the self-similar solutions. We note also 
that for C6 the k~ 3 ^ 2 scaling seems to end at the mag- 
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FIG. 6: Scaling as a function of Reynolds number Re of the 
maximum value of the dissipation over time, {top), of the time 
at which this maximum is reached (middle) and of the Taylor 
Reynolds number R\ computed at the instantaneous peak of 
dissipation (bottom); the straight line indicates the turbulent 
scaling R.\ ~ Re 1 / 2 . Symbols as in Fig. [5] 



netic Taylor scale, beyond which bending of magnetic 
field lines is felt and a steeper power law is possible, as 
already observed for a general non-symmetric flow [14| . 
but no such double power law is observed for the other 
two classes of flow. Furthermore, a bottleneck at the 
beginning of the dissipation range is noticebly absent or 
undetectable, likely due to the intrinsic nonlocality of 
nonlinear interactions in MHD (ifjj . 

Table |TT] gives quantitative data of the characteristic 
scales in the highest resolution runs of each class, includ- 
ing the integral and Taylor scales, as well as the Tay- 
lor Reynolds number, based on A T . When based on the 




FIG. 7: Total energy spectra averaged over At = 0.5 (1.5 
to 2 turnover times) around the time of maximum dissipation 
for different Reynolds numbers for the following flows (see Ta- 
ble U}: I runs compensated by k 2 (top), A runs compensated 
by fc 5//3 (middle) and C runs compensated by k 3 ^ 2 (bottom). 
Dash-triple dots, dash-dots, dashes, dots, and solid lines rep- 
resent respectively the runs 12 to 16, A2 to A6 and C2 to C6. 
Equivalent resolutions range from 128 3 to 2048 3 grid points. 

magnetic Taylor scale X M instead, the Taylor Reynolds 
numbers are 1600 (for 16), 950 (for A6), and 910 (for C6) 
respectively. 



V. DISCUSSION AND CONCLUSION 

There is a wealth of theoretical, phenomenological, ob- 
servational, and numerical studies of energy spectra for 
MHD turbulence. Solar wind data has seemed for a long 
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FIG. 8: ( Top) Magnetic to kinetic energy ratio as a function 
of wavenumber at later time (t £ [6, 6.5]) than in Fig. [3] and 
averaged over the same length of time. Again, we plot 16 
(solid), A6 (dashed) and C6 (dotted). Note that the excess 
of magnetic energy at low wavenumber in 16 has subsided to 
1/15 of its former value. {Bottom) Total energy spectra for the 
same runs, over the same time interval. The arrows indicate 
the new magnetic Taylor wavenumber for each run. In the 
insert, the same spectra are compensated by k 3 ^ 2 . Note the 
similar energy ratios and inertial range scaling for the three 
flows. 



while to favor a K41 classical spectrum, but a puzzling 
recent result is that, in some cases, both IK- and K41- 
type power laws have been observed for the velocity and 
magnetic field [Hj]; moreover, recent data on MHD tur- 
bulence in the plasma sheet using the CLUSTER suite 
of satellites seem to indicate that the inertial index of 
this flow varies but with more likelihood for a —2 law 
and to a lesser extent a —1.6 power law for the energy 
spectrum pfjj . The tendency toward K41 or IK dynam- 
ics has also been observed recently in several numerical 
simulations [HI, [HI, [23| with different forcing functions 
(see also [13 )• In fact, it has been shown numerically for 
the reduced MHD equations || that a whole palette of 
spectra is possible 0, [24| . 

While important theoretically for establishing how dy- 
namical energetic exchanges in MHD turbulence occur, 
distinguishing amongst possible spectral indices numeri- 
cally is difficult at best. Equally unclear is whether one 
power-law solution exists in all cases or several classes of 
universality are possible (48[, or if the resulting spectra 



can exhibit an arbitrary power law in a range constrained 
by integrability both at small and large scales, avoiding 
the so-called ultraviolet and infrared catastrophic diver- 
gences (49|. It was found recently that for a given set 
of initial conditions, an isotropic IK spectrum results in 
the large scales, with a transition at the magnetic Tay- 
lor scale towards a WT anisotropic spectrum, as men- 
tioned earlier, with both ranges only moderately resolved 
[l4| . Furthermore, boundary conditions and other forc- 
ing functions (e.g., random) may play a role as well. For 
example, it was found in a pioneering paper that the spec- 
trum (in reduced MHD) can become very steep when the 
time scale of the external forcing is varied [23 1 . 

To this debate we contribute here our numerical confir- 
mation that different energy spectra may indeed emerge 
in the specific case of decaying MHD turbulence in the 
absence of a uniform magnetic field. For this purpose, 
we have generalized the Taylor- Green flow to MHD and 
studied three different configurations that lead to three 
different types of behavior, even though, from a statis- 
tical point of view, these configurations were a priori 
equivalent since they had the same invariants (E T , H c 
and H M ) and the same equipartition between kinetic and 
magnetic energy at initial time. By taking advantage of 
the natural symmetries of these flows, we have been able 
to examine higher Reynolds numbers than for a full DNS 
for a given cost; we could also perform a convergence 
study in terms of Reynolds number, with R\ ~ 960 or 
above in all three flows at the highest resolution. We 
found that the I flow behaves somewhat differently, with 
a slower dissipation of energy at a given Reynolds num- 
ber and a lower maximum of dissipation, and that all 
three flows had different partition between the kinetic 
and magnetic energy at large scales, a partition that is 
likely at the origin of their different spectral behavior. 

That the observed discrepancy in spectral behavior is 
not linked intrinsically to a given type of initial condition 
but rather correlated with the ratio for time scales as 
displayed in Fig. U is confirmed by the following analysis. 
It is evident in Fig. [5] (bottom) that at late time (t ft 5.5) 
the ratio of magnetic to kinetic energy does not differ as 
much for the three types of flows studied in this paper. 
Therefore, we examine in Fig. [8] the ratio of magnetic 
to kinetic energy spectra (top) and the energy spectra 
(bottom) averaged over At = 0.5, as before, but now 
at a later time beginning at t = 6. Indeed, when the 
kinetic and magnetic energy are comparable, the energy 
spectra of the I and A flows likewise have comparable 
spectral index, the C flow (dotted line) being somewhat 
shallower although at this late time and lower Reynolds 
number, the difference is hard to tell. However, the I flow 
is clearly not as steep as at earlier times. 

The different power laws observed in this study can 
in principle be associated with known "solutions" (K41, 
IK, or WT) of MHD turbulence (omitting here intermit- 
tency corrections), and they are found to be correlated 
with the ratio of the nonlinear to the Alfvcn time. This 
is linked to the competition between nonlinear steepen- 
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ing and dispersion due to waves, which can interact as 
they propagate forwards and backwards along a mean 
field — here a local mean field, since there is no imposed 
background field. Whether the attractors for MHD (as- 
sociated with different spectral indices) are one, many, 
or infinite, remains to be determined. It is conceivable 
that multiple fixed points can co-exist, linked with K41 
(fluid- like) , IK (balance between steepening and disper- 
sion) and WT (turbulence and waves) dynamics. For 
example, one numerical study, though performed at low 
resolution, showed that there are indeed several possible 
attractive solutions for decaying incompressible MHD in 
the absence of an imposed magnetic field: one dominated 
by the velocity, another dominated by the magnetic field, 
and the third wit h q uasi-equipartition between the two 
modes of energy [50(. Are the solutions we observe in 
our study, which build up on a fast time scale of the or- 
der of the eddy-turnover time, associated with the slow 
dynamics of the attractive solutions of decaying MHD? 
This hypothesis might be tested by performing multiple 
selective decay simulations for very long times for the 
three TG MHD flows analyzed in this paper and variants 
thereof, a task left for future work. 

The initial conditions we have chosen are concen- 
trated in the largest scales in order to obtain the largest 
Reynolds numbers possible for these flows. Moving the 
characteristic wavenumber of the initial conditions to 
smaller scales would allow for an inverse transfer of mag- 
netic helicity to occur, and this in turn may well influ- 
ence the overall dynamics. A substantial amount of cross 
helicity, a different initial ratio of magnetic to kinetic en- 
ergy, or a non-unit magnetic Prandtl number could also 
have an important ciimpact on the resulting energy spec- 
tra. For example, it has been known for a long time that 
global and local correlation of velocity and magnetic field 
fluctuations can directly influence the spectral indices of 
the Elsasser variables and hence of all other spectra, with 
a steeper slope for the dominant mode. These effects have 
been studied in the context of isotropic turbulence clo- 
sures [5lj ; in weak turbulence [13]; in two-dimensional 
simulations [4l|, [HJ ; and in analytical studies of small- 
scale structures and wave interactions 0, HH . Similarly, 
it would be of interest to study the forced case, using the 
three (I, A, C) configurations studied in this paper. 

That the inertia! index of energy spectra can vary over 
a discrete or continuous range of values therefore has been 
supported by many previous studies. This lack of uni- 
versality has been studied analytically in terms of eddy 
shape and alignment with magnetic fluctuations (9j, nu- 
merically in terms of time scales in reduced MHD with 
forcing [231 ] . and experimentally in the different context 
of surface gravity wave turbulence [54| ■ It has also been 
linked to the amplitude of forcing in another experimen- 
tal study [55[ . Although the regime we study here numer- 
ically is one at high Reynolds number, where the role of 
waves is blurred by their interactions with eddies and also 
by their nonlinear steepening within structures, we nev- 
ertheless consider it an open problem to know whether 



other spectra than the three presented here can be found 
in MHD turbulence. Other open problems are to under- 
stand the interplay of eddies and waves in turbulence, 
particularly in the MHD case, and to unravel the role 
played by anisotropy induced by a large-scale magnetic 
field. 



Whatever power law may exist for an energy spectrum, 
it may moreover be affected additionally by the amount 
of intermittency present in the flow. Measurement of in- 
termittency (as performed recently for active regions of 
the sun |56|) in our flows — and comparison among the 
three classes — is left for future work. It should be linked 
to a study of the structures that develop, which both in- 
fluence and are influenced by the nonlinear dynamics. It 
is already known that current and vorticity sheets fold 
and roll-up at high Reynolds numbers [3, E3] ■ One fur- 
ther question is whether the set of anomalous exponents 
for the velocity and magnetic fields, or equivalently in 
terms of Elsasser variables, which differ at second order 
corresponding to the energy spectrum, have a common 
asymptote at high order, probably determined by the 
structures with steepest gradients. 



Further work should also include exploring runs at 
higher Reynolds numbers. Short of waiting for the next 
generation of resources (which will be made available 
to be developed, e.g. through the petascale computing 
initiative, one can resort to modeling methods, in ad- 
dition to direct numerical simulation. Akin to the one 
presented here, insofar as implementing a reduction of 
modes at a given Reynolds number, is the numerical 
algorithm that decimates modes (somewhat arbitrarily) 
in the dissipation range [57| . Another possibility is the 
use of large-eddy simulations that compare well against 
hig h-resolution direct numerical simulations, such as in 
[58l | . Of a different nature is the Lagrangian averaging 
approach, or alpha model, which can be viewed as a sort 
of direct numerical simulation methodology imposing a 
filter to the small scales by means of a closure consis- 
tent with preserving the Hamiltonian nature of the flow, 
although these averaged equations conserve the ideal in- 
variants using a different norm than £2 [59^ . With alpha 
models (see e.g. I60I) . it has recently been shown that 
the result found in of two inertial ranges for MHD, of 
isotropic IK followed by a weak turbulence spectrum, can 
be recovered at substantially lower cost. Using a combi- 
nation of such modeling tools may allow for parametric 
investigations of MHD turbulence and thereby lead to 
a better understanding of such flows as they occur in 
geospace, the hcliospherc, and the interstellar medium, 
and their influence for example on cosmic ray propaga- 
tion or on the solar-terrestrial interactions. 
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